Termination Mechanisms of Turing Patterns in Growing 
Systems 


Gabriel Morgado, Laurence Signon, Bogdan Nowakowski, Annie Lemarchand 


» To cite this version: 


Gabriel Morgado, Laurence Signon, Bogdan Nowakowski, Annie Lemarchand. Termination Mech- 
anisms of Turing Patterns in Growing Systems. Acta Physica Polonica B, Jagellonian University, 
Cracow, 2019, 50 (7), pp.1369. 10.5506/APhysPolB.50.1369 . hal-02324016 


HAL Id: hal-02324016 
https://hal.archives-ouvertes.fr/hal-02324016 
Submitted on 21 Oct 2019 


HAL is a multi-disciplinary open access 
archive for the deposit and dissemination of sci- 
entific research documents, whether they are pub- 
lished or not. The documents may come from 
teaching and research institutions in France or 
abroad, or from public or private research centers. 


L’archive ouverte pluridisciplinaire HAL, est 
destinée au dépôt et à la diffusion de documents 
scientifiques de niveau recherche, publiés ou non, 
émanant des établissements d’enseignement et de 
recherche français ou étrangers, des laboratoires 
publics ou privés. 


Termination mechanisms of ‘Turing patterns in 
growing systems 


Gabriel Morgado!, Laurence Signon?, Bogdan Nowakowski!?, and Annie 
Lemarchand? 


Institute of Physical Chemistry, Polish Academy of Sciences, Kasprzaka 
44/52, 01-224 Warsaw, Poland 
Sorbonne Université, Centre National de la Recherche Scientifique 
(CNRS), Laboratoire de Physique Théorique de la Matiére Condensée 
(LPTMC), 4 place Jussieu, case courrier 121, 75252 Paris CEDEX 05, 
France 
SWarsaw University of Life Sciences (SGGW), Department of Physics, 
02-776 Warsaw, Poland 


April 8, 2019 


Corresponding author: Annie Lemarchand, E-mail: anleQ@lptmc.jussieu.fr 


Abstract 


The question of the termination of a periodic spatial structure of Turing type ina 
growing system is addressed in a chemical engineering perspective and a biomimetic 
approach. The effects of the dynamical parameters on the stability and the wave- 
length of the structure are analytically studied and used to propose experimental 
conditions for which a Turing pattern stops by itself with a decreasing wavelength. 
The proposed mechanism is successfully checked by the numerical integration of the 
equations governing the dynamics of the activator and the inhibitor. We conclude 
that a local increase of the concentration of the reservoir which controls the injec- 
tion rate of the inhibitor into the system can be used to achieve the appropriate 
termination of a Turing pattern. 


1 Introduction 


During embryonic development, segmented structures of the body such as the spine and 
the digits are formed by the production of repeated elements. Since the seminal work 
of Turing [I| accounting for the formation of biological pattern in the framework of 
reaction-diffusion models, experimental evidences of Turing structures have been given 
in chemistry [2] [3] 4] and biology [5] [6]. Recent years have witnessed a growing num- 
ber of papers where reaction-diffusion principles are proposed to model the formation 
of biological periodic spatial structures [7 [8] [9] ZO, (11) 12) [13]. Following Turing, a 
two-component chemical system composed of an autocatalytically-produced activator by 
consumption of an inhibitor that diffuses faster may produce periodic spatial structures 
such as stripes in one-dimensional (1D) systems and hexagons in 2D. In other words, a 
Turing pattern emerges by local self-activation and lateral inhibition [14]. The concepts 
developed to model living systems provide a source of inspiration in chemical engineering 
[15] 16} [A C8, [19] 21, LA. However, standard models of Turing patterns generate struc- 
tures in infinite systems and the question of the termination of a striped structure in a 
finite system arises in a perspective of biomimetism in material science. Specifically, it 
is desirable to find experimentally achievable conditions creating a finite-size structure, 
whose growth stops by itself with decreasing oscillation amplitude and respects the de- 
crease of the wavelength during the termination process. To this aim, we use an as simple 
as possible reaction-diffusion model [23] admitting a Turing structure developing behind 
a propagating wave front and examine the effect of all parameters on both the stability 
and the wavelength of the structure [5] [22]. We already used a stochastic approach to a 
Turing pattern [23] and showed that, contrary to intuition, the internal fluctuations may 
have a constructive effect and stabilize the structure outside the domain of stability pre- 
dicted by a deterministic description. Here, we adopt a macroscopic approach. Our goal 
is to select suitable conditions from this systematic approach and to propose termination 


mechanisms compatible with processing in chemical engineering. 


The paper is organized as follows. In section 2, a reaction-diffusion model involving 


local activation and long-range inhibition is presented. An analytical stability condition 
and the wavelength expression of the Turing pattern are given. The influence of the 
parameters of the model on the stability and the wavelength of the pattern are studied in 
section 3. The analysis of the results leads to the selection of a user-friendly termination 
mechanism in the framework of chemical engineering. The analytical predictions regarding 
stability and wavelength are compared to numerical results for the chosen mechanism. 
Section 4 contains conclusions. The possibility that the different mechanisms exhibited 


could be found as termination scenarios in morphogenesis is raised. 


2 Model 


We consider the following reaction mechanism inspired from the Schnakenberg model [24] 


and the Gray-Scott model [25] 


2A +B — 3A (2) 


k3 
k_3 


where R; and Rs are reservoirs. The concentrations À; and Ry of the reservoirs are 
assumed to be constant in time. The reaction given in Eq. (2) autocatalytically produces 
species A and consumes species B. Due to the ability of accelerating its own production, 
species A is called an activator whereas species B, which is consumed by the same process, 
is called an inhibitor. The macroscopic dynamics of the system is governed by two partial 


differential equations [9] [23] 


OA OA 
OE = —kA+kA?B + Daas (4) 
ðB B 


ot = k_3Ro = k3B = kə A? B + Dega (5) 


for the concentrations A and B of the activator and the inhibitor supposed to have different 


diffusion coefficients D4 and Dg. For appropriate rate constant values, such that 
A = (k_3R2)° = Ak? k3/ke > 0 (6) 
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Figure 1: Largest eigenvalue u of the linear operator M versus square of Fourier mode q?. 
Solid line: k_3R2 = 8.76. Dashed line: k_3R2 = 10. Other parameter values: kı = 2.92, 
ke = 1, hy = 2.19, Dg =1, Dp = 10. 


the system admits two steady states (Ap = 0, Bo = k_3R2/k3) and 


_ k3R+VA 
Ar = >) — (7) 
Bp = k_3Rə — VA (8) 
2k3 


that are stable with respect to homogeneous perturbations. A linear stability analysis of 
Eqs. (45) reveals that the steady state (Ar, Br) can be destabilized by inhomogeneous 
perturbations 5,810123]. The Fourier transforms A,(t) = [°° A(x,t)e “dx and B,(t) = 
Jo, B(x, the dz, where q is the Fourier mode, are introduced. In the Fourier space, the 


linear stability operator M is given by: 


o kı — Dag? k> AF 
M = | —2k — 28 47 = Dee? (9) 


The eigenvalues of the matrix M obey the characteristic equation u? + au + 8 = 0, with 
a = kı — P Ar- (Da + Dp) and 8 = 2k? Ar/Br — (kı — Daq?)(k_-3R2/Br + Dpq’). 
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The Turing structure develops if the largest eigenvalue 


1 | kok_3Ro 
E Cr 
1 


5 Ar — (Da + Dg) + 


kok_3R | 
(r + Ap + (De — Dan) — Sky ky AB, (10) 
al 


is real and positive [5] B]. Indeed, a system of differential equations, linearized around 
a homogeneous steady state, is easily solved by diagonalizing the linear operator. Then, 
the solution is a linear combination of eigenmodes which exponentially depend on time 
according to the corresponding eigenvalues. A term associated with a real, positive eigen- 
value grows in time, leading to trajectories in the concentration space that move away 
from the fixed point [5]. In the studied system, the destabilization of the steady state 
occurs in favor of a Turing pattern. Equation (10) imposes conditions on the parameter 
values. In particular, the diffusion coefficient Dg of the inhibitor B must be sufficiently 
larger than the diffusion coefficient D4 of the activator A: The destabilization of the homo- 
geneous steady state (Ar, Br) requires local self-activation, ensured by the autocatalytic 
production of the activator through the reaction given in Eq. (2), as well as long-range 
inhibition, due to the larger diffusion coefficient of the inhibitor. The mode quax, which 


maximizes the eigenvalue y, is the most unstable Fourier mode: 


B Ar(Da + Dg)y/2kık2D4/Dpg — kı — kək-3R2Ar/kı (11) 


— Ds- Da 
In order to account for the termination of the Turing pattern in a growing system, 
including the fact that the structure ends with a gradually shorter spatial oscillation, 
we need to find conditions for which the structure tends to lose its stability while its 


wavelength decreases. The wavelength of the periodic structure is given by: 


27 


= (12) 


max 


Turing structure becomes unstable as the value of largest eigenvalue vanishes for the mode 


max associated with the maximum of ju: 


Umax = =F Ar — (Da + Dp) @raxt 
2 ky 
Lis ? 
(r + LA + (Dg -— Das) — 8kik2 Ap (13) 
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Figure 2: Top: Scaled wavelength \/Az of Turing pattern versus diffusion coefficient D 4 
of species A. Bottom: Maximum value [imax of the largest eigenvalue of the linear operator 
M versus D1. Parameter values: kı = 2.92, kə = 1, k3 = 2.19, k_3R2 = 8.76, Dg = 10, 
Arx = 0.1. 
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Figure 3: Top: Scaled wavelength A/Az of Turing pattern versus diffusion coefficient Dz 
of species B. Bottom: Maximum value [max of the largest eigenvalue of the linear operator 
M versus Dg. Parameter values: kj = 2.92, kg = 1, k3 = 2.19, k_3R2 = 8.76, Da = 1, 
Ax = 0.1, 


with dmax given in Eq. (11). Figure 1 illustrates the behavior of Hmax for parameter values 
associated with a stable Turing pattern with Umax > 0. It is also shown that it is sufficient 
to increase the value of k_3R, to shift the curve u(q?) down and lose the stability of the 


Turing pattern. 


In the next section we investigate the behavior of À and [max as each parameter 
controlling dynamics varies. Specifically, we aim at identifying diffusion coefficients or rate 
constants whose variation leads both to a decrease of the wavelength and a destabilization 


of the Turing structure, i.e. negative values for the maximum of the eigenvalue. 


3 Results 


The concentration Rə of the inhibitor reservoir is first assumed to be homogeneous. Fig- 


ures 2 and 3 show the variations of the wavelength À and the maximum value Hmax of the 
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eigenvalue with one of the diffusion coefficients, the other parameters being constant. The 
results are deduced from Eqs. (IT) and (12) for À and Eq. (18) for Umax, the expressions 
of the steady state (Ar, Br) being given in Eqs. (7) and (8). To facilitate the comparison 
with the numerical integration of Eqs. (4) and () that will be performed in the following, 
the wavelength is given in number of spatial cells of length Ag = 0.1. As shown in Fig. 
2, the decrease in the maximum value [max of the eigenvalue as the diffusion coefficient 
D4 of species A increases is accompanied by an increase of the wavelength À: The loss of 
stability of the Turing structure occurs with an increase of the spatial period. We con- 
clude that a variation of the diffusion coefficient D4 cannot be argued as a justification 
of the termination process. The behavior with respect to the diffusion coefficient Dp of 
species B is different. The simultaneous loss of stability of the structure and the decrease 
of the wavelength are observed in Fig. 3 as Dg decreases: The diffusion coefficient Dp of 


species B can be considered as a suitable parameter in the search for a termination model. 


Figures 4-7 show the variations of the wavelength À and the maximum value [max of 
the eigenvalue with rate constants. The variations of À are given in a bounded interval 
of rate constant values, in which the Turing pattern is stable. At one of the endpoints 
of the interval, the eigenvalue Hmax vanishes and at the other endpoint, the condition of 
existence of the steady state (Ar, Br) given in Eq. (6) is no longer satisfied. The two 
desired behaviors, i.e. the decrease of both À and max, are observed as kı decreases, kə 
increases, k decreases, and k_3R, increases. For an assumed homogeneous concentration 
Rə of the reservoir, the variations of À and [imax with Rə are analogous to the variations 
with k_3R2. According to the chemical reaction given in Eq. (1), decreasing the rate con- 
stant kı tends to increase the concentration of species A. Following Eq. (2), increasing the 
rate constant kz of the autocatalytic step tends to increase the concentration of species A 
and decrease the concentration of species B. This last result seems to be inconsistent with 
the consequences drawn from the decrease in k or the increase in k_3R2, which result 
in increasing the concentration of species B according to Eq. (B). However, we already 


stated that increasing B through soliciting the reservoir Rə results in consuming species 
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Figure 4: Top: Scaled wavelength \/Azx of Turing pattern versus 


rate constant kı of 


the chemical reaction given in Eq. (). Bottom: Maximum value pmax of the largest 


eigenvalue of the linear operator M versus kı. Parameter values: 
k_3Ro = 8.76, Da = 1, Dz = 10, Ar = 0.1. 
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Figure 5: Top: Scaled wavelength \/Az of Turing pattern versus rate constant k2 of 
the chemical reaction given in Eq. (2). Bottom: Maximum value [imax of the largest 
eigenvalue of the linear operator M versus ky. Parameter values: kı = 2.92, k3 = 2.19, 
k_3Ro = 8.76, DA = 1, Dp = 10, Ag= 0.1, 


B faster by the autocatalytic step given in Eq. (2) [9) 26]. In particular, we observed 
that introducing a local source of species B leads to the nonintuitive local decrease of B 
concentration. Hence, all the variations of the rate constants that lead to a loss of stability 
of the Turing pattern are eventually associated with an increase of A concentration and 


a decrease of B concentration. 


The diffusion coefficients and the rate constants characterize dynamics and are intrin- 
sic to the reaction-diffusion system. Nevertheless, it is always possible to imagine spatial 
variations of the dynamical parameters. Well-chosen variations of the diffusion coefficient 
D; of the inhibitor and each of the four rate constants of the chemical mechanism could 
be a priori used to build a termination model. In the framework of the application to 
developmental biology, steric hindrance and molecular crowding in the tail of an embryo 


may be invoked to justify the decrease of the diffusion coefficients. In chemical engineer- 
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Figure 6: Top: Scaled wavelength A/Az of Turing pattern versus rate constant k3 of the 
forward chemical reaction given in Eq. (8). Bottom: Maximum value pmax of the largest 
eigenvalue of the linear operator M versus k3. Parameter values: kı = 2.92, kə = 1, 
k-3Rə = 8.76, Da = 1, Dp = 10, Ar = 0.1. 


ing, a local increase of temperature could be used to induce a local increase of the rate 
constants. However, local increase of confinement or temperature is susceptible to simul- 
taneously affect several dynamical parameters [13] 22, 27, [28] 29}, [30] B1 [82]. Whereas 
a decrease of Dg is desired to destabilize the Turing pattern while decreasing its wave- 
length, a simultaneous decrease of D4 would be detrimental. Similarly, an increase of ko 
and k_3 due to temperature increase could be satisfying but the joint decrease of kı and 
k3 could blur the effect on the Turing structure. The simplest way to imagine the control 
of a targeted parameter leading to the desired behavior is to impose well-chosen spatial 
variations of the reservoir concentration Rə. Indeed, the product k_3R2 plays the role of 
an apparent rate constant for the backward reaction given in Eq. (B) that can be fixed 
at will in chemical engineering in the case of a single dynamical system with uniquely 


defined intrinsic parameters. 
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According to Fig. 7, increasing Rə tends to destabilize the Turing pattern and decrease 
its wavelength. We examine if the results deduced from a stability analysis can be used 
in a dynamical approach. The results of the numerical integration of Eqs. (4) and (6) for 
a homogeneous concentration Rə and a piecewise linear profile are given in Fig. 8. The 
initial condition is a step function between the steady state (Ar, Br) in the first 10 cells 
on the left and the steady state (Ao, Bo) in the remaining cells. The initial total number 
of cells is set at no = 610. Introducing the cell label i = x/Ax, where Az is the cell 


length, and the discrete time s = t/At, where At is the time step, we choose: 


A(i, s = 0) = Ar, B(i,s = 0) = Br, for 1<4< 10 (14) 


A(i, s = 0) = Ap, B(s,s = 0) = Bo, for 11 <i < no (15) 


To account for the growth of the system and simultaneously avoid boundary effects that 
may alter the wavelength of the structure [16], we impose a fixed boundary on the left 
and a free growing end on the right [9] 23, 26]. For parameter values for which the steady 
state (Ar, Br) is unstable with respect to inhomogeneous perturbations, a Turing pattern 
develops after the passage of a chemical wave front. More precisely, according to Eqs. 
(4) and (5) and due to the no-flux boundary conditions applied on the left boundary, the 


concentrations in the first cell obey 


A(1,s+1) = A(1,s) — kAtA(1,s) + k2AtA(1, s} B(1, s) + 
Darf lAl, 8) ~ AGL) (16) 
Bil,s+1) = B(1,s)+k_3R,At — k3AtB(1, s) — k,AtA(1, s)*B(1,s) + 
At 
Darg BO, s) — B(1,s)) (17) 


so that both A(1,s) and B(1, s) are extremum of the Turing pattern in the first spatial 


cell 7 = 1. 


Spatial cells are added to the right end of the system at the front speed to counter- 
balance the progression of the wave front and mimic system growth: At all the discrete 


times s for which the concentration B(n — 600, s) of species B in the n — 600 cell becomes 
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Figure 7: Top: Scaled wavelength À/Ax of Turing pattern versus rate constant k_3R2 of 
the backward chemical reaction given in Eq. (8). Bottom: Maximum value [max of the 
largest eigenvalue of the linear operator M versus k_3R2. Parameter values: kı = 2.92, 
ko = 1, kg = 2.19, Da = 1, Dp = 10, Ar = 0.1. 


smaller than 0.99B,, the total number n of cells is increased by 1. Provided that the 
front propagates at a speed smaller than Ax/At, this protocol ensures that a layer of 
about 600 cells remains in the stationary state (Ag, Bo) on the right of the system, so that 
the propagation of the front is not significantly affected by the finite size of the system. 
To draw Fig. 8b, we have chosen the parameter values given in the caption of Figs. 1 
and imposed k_3 = 8.76 for the following spatial profile for the concentration Ry of the 


inhibitor reservoir: 


Rə = 1,  forl<i< 500 (18) 
Ry = 2.83 x 10 %+0.858,  for500 < i< 1000 (19) 
Rə = 1.14,  for1000 <i (20) 


The simulation is stopped at time tena for which the wave front has passed cell à = 1000. 
It is worth noting that the Turing pattern is unchanged for larger values of the final inte- 


gration time. Then, only the position of the concentration gradients associated with the 
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traveling wave evolve in time but the Turing structure has stopped growing and remains 
in a steady state with a fixed number of wavelengths. As desired, the increase of the 


concentration Rg leads to the termination of the Turing structure. 


As illustrated in Fig. 7, the Turing structure is expected to be stable in the range 
1 <7 < 500 for which k-3R2 = 8.76 and unstable in the range i > 1000 for which 
k_3R2 = 10. More precisely, according to Eq. (18), the maximum of the eigenvalue pmax 
vanishes for k_3R2 = 9.75, i.e. Rə = 1.11 for k_3 = 8.76, which occurs in spatial cell 
i = 900. Hence, the Turing pattern is predicted to be stable in the range 0 < à < 900 
and unstable beyond this domain. The results shown in Fig. 8b confirm the analytical 
predictions. The amplitude of the spatial oscillations decreases between à © 500 and 
i ~ 1000. The system is in a steady state in the range 1000 < à < 1500. 
The increase of Ra not only destabilizes the Turing structure but also modifies the steady 
state values and the propagation speed of the wave front. The comparison between Figs. 
8a and 8b shows that, as Rə increases, the wave front propagates faster, Ar increases, 
Br decreases and Bp increases. As a consequence of the variation of Ar and Br, the 
oscillations of A and B concentrations are not symmetrical in the range 500 < i < 
900. The decrease of the wavelength predicted in Fig. 7 is more difficult to check by 
a qualitative analysis. Using the numerical results illustrated in Fig. 8b, we evaluate 
the local wavelength by computing the number of cells between two minima of the A 
concentration profile. The results are given in Fig. 9 and compared to the analytical 
prediction deduced from Eqs. (LJ) and (12). The agreement between the numerical and 
analytical results is very satisfying in the range 600 < 4 < 900. Oscillations of very small 
amplitude are observed in Fig. 8b in the range 900 < à < 1000, proving that a very 
damped Turing structure remains in a small area where instability was predicted. The 
wavelength of the structure in the range 1 < à < 500 is slightly affected by the increase 
of Ry from cell i = 500 but the deviation from the analytical prediction is only 2.5 
percent. This small difference is related to the linear approximation used in wavelength 


evaluation that neglects non-linear terms that may be more important for large structures. 
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Interestingly, the wavelength is sensitively decreased in the expected area in which the 
concentration of the reservoir Ry increases: As shown in Fig. 9, the wavelength is reduced 
from 72 spatial cells to less than 61, before the structure disappears. We conclude that 
an increase in the concentration of the reservoir Rə related to the inhibitor B is sufficient 
to account for the destabilization of the Turing pattern associated with a decrease of the 
wavelength. As anticipated by the results given in Fig. 7, according to which an increase 
of Rə decreases the wavelength À and leads to a negative eigenvalue [max around (Ar, Br), 
we suggest that an appropriate spatial variation of Rə can be used in chemical engineering 
to stabilize the homogeneous steady state and induce a termination of the Turing pattern 


in a growing system. 


4 Conclusion 


In a biomimetic approach, we have addressed the question of the termination of a Tur- 
ing structure in a growing system. A free boundary is imposed at the growing part, 
which ensures that the wavelength of the pattern is not perturbed by fixed boundary 
conditions. After deriving analytical expressions for the stability condition and the wave- 
length of the structure, we perform a systematic analysis of the effect of all dynamical 
parameters on the pattern. Apart from the variation of the diffusion coefficient of the ac- 
tivator, a well-chosen variation of the dynamical parameters leads to the desired behavior, 
i.e. the simultaneous loss of stability and the decrease of the wavelength. In particular, 
an increase of the effective rate constant k_3R2, where k_3 is the rate constant of the 
reaction injecting the inhibitor from the reservoir at the concentration Rə, is associated 


with a destabilization of the Turing pattern accompanied by a decrease of the wavelength. 


Imposing a spatial variation of the concentration of the reservoir Rə turns out to be an 
appropriate protocol for chemical engineering. However, the proposed procedure imposes 
the total length of the structure but not its number of wavelengths. In the framework of 
developmental biology, for example in the case of the growth of the digits or the spine 


of the vertebrates, the termination process has to respect the total number of segments 
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Figure 8: Spatial profiles deduced from the numerical integration of Eqs. (4) and (5) for 
kı = 2.92, kə = 1, ka = 2.19, he = 8.76, Da = 1, Dg = 10, At = 10, tena/ At = 
2000000, Az = 0.1. Black dotted line: Imposed concentration Rə of the reservoir. (a): 
Homogeneous concentration Rə = 1, (b): Piecewise linear Rọ profile. Blue dashed line: 
Concentration of species A versus cell label x/Ax. Red solid line: Concentration of species 
B versus cell label x/Ax. 
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Figure 9: Spatial variation of the scaled wavelength A/Azx of Turing pattern versus cell 
label x/Ax for the parameter values given in the caption of Fig. 8b. Symbols: Results 
deduced from the numerical integration of Eqs. (4) and (5). Solid line: Analytical 
prediction given in Eq. (12). 
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for a possible variation in the length of the global structure. Therefore, it is necessary to 
imagine that the system itself is able to count the number of already formed segments and 
to trigger the variation of a parameter leading to smaller subsequently formed segments. 
If the concept of the Turing structure is kept in the formation of biological patterns, the 
presented results could be used to suggest such relevant parameters. The local increase 
of the rate constant k_3 that would be activated when a given number of segments has 
already been formed can be straightforwardly proposed. Similarly, the local increase of 
the rate constant kə controlling the autocatalytic production of the activator or the local 
decrease of the rate constant kı or k3, associated with the absorption of the activator or 
the inhibitor by reservoirs, would lead to the desired phenomenon. The local decrease of 
the diffusion coefficient of the inhibitor offers an alternative. The nature of the mechanism 
that would trigger such a response of the system when a given number of segments has 


been created remains an open question. 
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